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Abstract. We consider the motion of point masses given by a natural extension 
of Newtonian gravitation to spaces of constant positive curvature. Our goal is 
to explore the spectral stability of tetrahedral orbits of the corresponding 4-body 
problem in the 2-dimensional case, a situation that can be reduced to studying 
the motion of the bodies on the unit sphere. We first perform some extensive 
and highly precise numerical experiments to find the likely regions of stability and 
instability, relative to the values of the masses and to the latitude of the position of 
three equal masses. Then we support the numerical evidence with rigorous analytic 
proofs in the vicinity of some limit cases in which certain masses are either very 
large or negligible, or the latitude is close to zero. 
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1. Introduction 

The goal of this paper is to study the spectral stability of tetrahedral orbits in 
the 2-dimensional positively curved 4-body problem, i.e. when four point particles 
of positive masses move on the unit sphere according to a gravitational law that 
naturally extends the Newtonian potential to spaces of constant curvature. This is 
a particular case of the curved n-body problem, n > 2, which has been only recently 
derived in a suitable setting, both for constant positive curvature (i.e. 2- and 3- 
dimensional spheres) and for constant negative curvature (i.e. 2- and 3-dimensional 
hyperbolic spheres), [12], [13], [6], [7]. 

The case n = 2 dates back to the 1830s, when Janos Bolyai and Nikolai Lobachevsky 
introduced it for the hyperbolic space H^, [2], [28]. The equations of motion for 
constant nonzero curvature are given by the cotangent of the spherical distance, in 
the positive case, and the hyperbolic cotangent of the hyperbolic distance when the 
space curvature is negative. For zero curvature, the classical Newtonian equations of 
the n-body problem are recovered. The analytic expression of the potential is due to 
Ernest Schering, [31], [32], for negative curvature, and to Wilhelm Killing, [17], [18], 
[19], for positive curvature. The problem also became established due to the results of 
Heinrich Liebmann, [23], [21], [25]. The attempts to extend the problem to spaces of 
variable curvature started with Tullio Levi-Civita, [21], [22], Albert Einstein, Leopold 
Infeld, Banesh Hoffman, [15], and Vladimir Fock, [16], and led to the equations of 
the post-Newtonian approximation, which are useful in many applications, including 
the global positioning system. But unlike in the case of constant curvature, these 
equations are too large and complicated to allow an analytic approach. 

It is important to ask why the above extension of the Newtonian potential to 
spaces of constant curvature is natural, since there is no unique way of generalizing 
the classical equations of motion in order to recover them when the space in which 
the bodies move flattens out. The reason is that the cotangent potential is, so far, 
the only one known to satisfy the same basic properties as the Newtonian potential in 
its simplest possible setting, that of one body moving around a fixed centre, the so- 
called Kepler problem [20] • Two basic properties stick out in this case: the potential 
of the classical Kepler problem is a harmonic function in M^, i.e. it satisfies Laplace's 
equation, and it generates a central field in which all bounded orbits are closed, a 
result proved by Joseph Louis Bertrand in 1873, [1]. 

On one hand, the cotangent potential approaches the classical Newtonian potential 
when the curvature tends to zero, whether through positive or negative values. On 
the other hand, the cotangent potential satisfies Bertrand's property for the curved 
Kepler problem and is a solution of the Laplace-Beltrami equation, [7], [20], the 
natural generalization of Laplace's equation to Riemannian and pseudo-Riemannian 
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manifolds, which include the spaces of constant positive curvature k > we are 
interested in here. 

In the Euclidean case, the Kepler problem and the 2-body problem are equivalent. 
The reason for this overlap is the existence of the linear momentum and centre of 
mass integrals. It can be shown with their help that the equations of motion are 
identical, whether the origin of the coordinate system is fixed at the centre of mass 
or fixed at one of the two bodies. For nonzero curvature, however, things change. 
The equations of motion of the curved n-body problem lack the linear momentum 
and centre of mass integrals, which prove to characterize only the Euclidean case, [7j, 
[8], [12]. Consequently the curved Kepler problem and the curved 2-body problem 
are not equivalent anymore. It turns out that, as in the Euclidean case, the curved 
Kepler problem is Liouville integrable, but, unlike in the Euclidean case, the curved 

2- body problem is not, [53], [3S], [SZ]- As expected, the curved n-body problem is 
not integrable for n > 3, a property also known to be true in the Euclidean case. 

A detailed bibliography and a history of these developments appear in [7] . Notice 
also that the study we perform here in §^ is not restrictive since the qualitative 
behaviour of the orbits is independent of the value of the positive curvature, [7|, [T2] . 

The current paper is a natural continuation of some ideas developed in [29], which 
studied the stability of Lagrangian orbits (rotating equilateral triangles) of the curved 

3- body problem on the unit sphere, S^, both when the mutual distances remain con- 
stant and when they vary in time. The former orbits, called relative equilibria, are a 
particular case of the latter, and they are part of the backbone towards understanding 
the equations of motion in the dynamics of particle systems, [7], [TU]. Unlike in the 
classical Newtonian 3-body problem, where the motion of Lagrangian orbits takes 
place in the Euclidean plane, the Lagrangian orbits of exist only when the three 
masses are equal, [12], [7]. But equal-mass classical Lagrangian orbits are known to 
be unstable, so it was quite a surprise to discover that, in S^, the Lagrangian relative 
equilibria exhibit two zones of linear stability. This does not seem to be the case 
for constant negative curvature, i.e. in the hyperbolic plane H^, as some preliminary 
numerical experiments show. Consequently, the shape of the physical space has a 
strong influence over particle dynamics, therefore studies in this direction promise 
to lead to new connections between the classical and the curved n-body problem. 

The result obtained in [22] thus opened the door to investigations into the stability 
of other orbits characteristic to S^, and tetrahedral solutions came as a first natural 
choice, since the experience accumulated in the previous study could be used in this 
direction, as we will actually do here. 

The paper is organized as follows. In Section 2, we introduce the tetrahedral 
solutions in S^, i.e. orbits of the 4-body problem with one body of mass mi fixed at 
the north pole and the other three bodies of equal mass m located at the vertices of 
a rotating equilateral triangle orthogonal to the 2;-axis. If the triangle is above the 
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equator, i.e. the z coordinate of the three equal masses is positive, the tetrahedral 
relative equilibria exist for any given masses. If the triangle is below the equator, the 
relative equilibria exist just for some values of the masses. To approach the spectral 
stability of the relative equilibria, we compute in Section 3 the Jacobian matrix of 
the vector field at the relative equilibria, which become fixed points in the rotating 
frame. 

The study of the stability starts in Section 4, where we analyze three limit prob- 
lems, first taking the mass at the north pole mi = 0, i.e. F := m\jm = 0. In this 
case, the tetrahedral relative equilibria are spectrally stable for z < and unstable 
for z > 0. In the second limit problem the mass mi is very large when compared 
to the other three masses. Taking e := l/T = m/mi, in the limit case e — j- 0, the 
problem reduces to three copies of 2-body problems, formed for the mass at the north 
pole and a body of zero mass. The changes in this degenerate situation for small 
e > are studied in Section 6, where we also consider the third limit problem, for 
which we take z = and let the parameters F or e and z move away from zero. To 
reach this point, we previously perform in Section 5 a deep and highly precise nu- 
merical analysis to determine the regions of stability according to the values of z and 
of the masses. Our main results occur in Section 6, where using the Newton polygon 
(including the degenerate cases) and the Implicit Function Theorem we study all 
bifurcations that appear when the limit problems are perturbed and draw rigorously 
proved conclusions about the spectral stability of tetrahedral relative equilibria. We 
end this paper with a full bifurcation diagram and an outline of future research 
perspectives. 



2. Tetrahedral orbits in 

Consider four bodies of masses mi,m2,m3,m4 > moving on the unit sphere S^, 
which has constant curvature 1. Then the natural extension of Newton's equations 
of motion from to is given by 

(1) q,= "^i^"^ 7 ^"^^ ' ^i'^^ - ■ ^M^. (q.-q.) = l, ^ = 1,2,3,4, 

where the vector = {xi,yi,Zi) gives the position of the body of mass mi,i = 
1, 2, 3, 4, and the dot, ■ , denotes the standard scalar product of M^, [13], [30]. These 
equations are known to be Hamiltonian, [7]. 

By a tetrahedral solution we mean an orbit in which one body, say mi, is fixed 
at the north pole (0,0, 1), while the other bodies, m2 = = m^^ =: m, lie at the 
vertices of an equilateral triangle that rotates uniformly in a plane parallel with the 
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equator z = 0. In other words, we are interested in solutions of the form 





= 0, 


yi 


= 0, 




Zl 


= 1, 






X2 


= r cos cut, 


2/2 


= r sin ut, 




Z2 


= ±( 




)l/2 


X3 


= r cos{ut + 27r/3), 


ys 


= r sin (cut - 


h27r/3), 


Z3 


= ±( 




J 1/2 


X4 


= r cos(cjt + 47r/3), 




= r sin(a;t - 


h47r/3), 


Zi 


= ±( 




1/2 



where r and u are constant, r denotes the radius of the circle in which the triangle 
rotates, and u represents the angular velocity of the rotation. A straightforward 
computation shows that 

2 24m mi 

^ ^ r3(12 - 9r2)3/2 ^3(i _ ^2)1/2 

where we take the plus or the minus sign depending on whether z := Z2 = = z^ is 
positive or negative, respectively. 

The purpose of this paper is to provide a complete study of the spectral stability of 
such orbits, which are obviously periodic. For this, we will use rotating coordinates, 
in which the above periodic relative equilibria become fixed points for the equations 
of motion. Recall that a fixed point is linearly stable if all orbits of the tangent flow 
are bounded for all time, and it is spectrally stable if no eigenvalue is positive or has 
positive real part. Linear stability imphes spectral stability, but not the other way 
around. Nevertheless, spectral stability fails to imply linear stability only in the case 
of matrices with multiple eigenvalues whose associated Jordan block is not diagonal. 

To achieve our goal, we further consider the coordinate and time-rescaling trans- 
formations 

= (Xj, Vi, Zi) Qi = {Xi, Yi), t = r^/'^T, 

x,=rX,, yi = rY,, Zi = ±[l - r^X^ + Y^)Y'\ i = 1,2,^ A- 
A simple computation shows that if we choose ut = Qt, the angular velocity relative 
to the new time variable r takes the form 

24m mi ^^"^ 

;i2 -9r2)3/2 ^ (l_r2)V2 

With the above transformations, and using the fact that 



n = ± 



X[ = r^l^x,, Yl = X'l = r%, F/' = r%, z = 1,2,3,4, 

the equations of motion become 



4 



Q:: = r' E "n^^ ,.{y?^^ -r^(q.-q.)Q. ^ = 1,2,3,4, 
where ' = ^ and /jj = (q^ ■ q^) =r'^{XiXj + YiYj) + ZiZj. 
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We further introduce the rotating coordinates ^j, rji, 2 = 1,2, 3, 4, with 



Yi 



T 



Then + r]ir]j = XiXj 
when i = j. Moreover, 



z = 1,2,3,4, where 7^(^]r) = 
YiYj, i,j e {1,2,3,4}, expressions that take the value 1 



cos Qt 
sin Qt 



sin Qt 
cos f2r 



2VL 



Vi 



Y!' 



1,2,3,4. 



A straightforward computation shows that the new equations of motion have the 
form 



(2) 



2Vt 



Vi 



where 

(3) 



24m 



± 



n9i 



Vj 



■k 



nil 



60 + ViVj 
9i,j 



^12-9r2)3/2 - (1 _ ^2)1/2' 
2 /■2 I „ 



(l-rVD(l-rV^) 



l\2 



1 — r'^p^ 



fi,j = r^i^i^j + ViVj) + ZiZj = r'^iiiij + y]ir]j) + s 



sign(2;i2j) 



sign(zj 



s, z = 1 or j = 1, 

which imphes that z = sa/1 — r"^- 

Before we start to study the stabihty of the tetrahedral relative equilibria, we 
must see for what values of the masses they exist. For this purpose, we will prove 
the following result. 

Proposition 1. Consider a tetrahedral orbit of the curved A-body problem in with 
the mass nii > fixed at the north pole (0, 0, 1) and the masses m2 = = 777,4 =: 
m > fixed at the vertices of an equilateral triangle that rotates uniformly on in a 
plane parallel with the equator z = 0. Then, if the triangle is above the equator, i.e. 
< z < 1, tetrahedral relative equilibria exist for any given masses. If the triangle 
is below the equator, i.e. —l<z<0, then 

(i) if < mi < for any positive value of up to a maximum it can attain, 
tetrahedral relative equilibria exist. In this case, if < mi < for any positive 
value of Vt^ , such that < < ^ — there is a unique z with — 1 < 2 < 
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corresponding to a relative equilibrium. Otherwise, there are two distinct values of 
z G (—1,0), each corresponding to a different relative equilibrium; 
(a) if nil > then there are no tetrahedral relative equilibria. 



Proof. From equation ([3]), a tetrahedral relative equilibrium must satisfy the condi- 
tion 

ft^ = m I — I ■= I , where z = ±Vl — r^, T = 

z 73(1 + 3^2)3/2; rn 




Figure 1. Plot of F(m, T),u e [0, 1], for six distinct values of F. From 
top to bottom, we took F = 0.3, 1/^3, 0.7, 16/(9^3), 1.2 and 1.5, 
respectively. 

For < z < 1, since the right hand side is positive, the statement in the proposition 
is obvious. For — 1 < 2 < 0, we introduce u = —z. Then the equation above can be 
written as 

(4) n^ = mF(u;F), 

where 

F[u- F) = -- + — ? -. 

^ ' U 73(1 + 3^2)3/2 

The behaviour of F[u] F) for different values of F is summarized in the Figure [1] The 
critical points of F{u\ F) satisfy F = g{u) where 

72^3 

9W 



\/3(l + 3^2)5/2 • 
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The function g{u) has a maximum at n = l/-\/2 equal to 144/(25-\/l5). For any 
< r < 144/(25\/T5) there exist values < ui < 1/v^ < U2 such that ^ = 0. It is 
clear that F has a maximum at Ui and a minimum at U2- Moreover, 

r(i - 6^^) 

and 

If r > ^ then Ml > 1/V6. Therefore F{ui,r) < and F{u,r) < for any 

u G (0, 1). Therefore equation (jlj) has no real solutions and there are no tetrahedral 
relative equilibria. If F < then F(ui,F) > and for any value of smaller 

than mF{ui;T) equation has real solutions. In this case, if < F < ^ and 

fi^ < F(1;F) = l/VS — F, there is a unique real solution u E (0,1), otherwise 
there are exactly two real solutions < uri < < ^ that correspond to relative 
equilibria. 

So, if I/VS < F < 16/(9^3), we can say that for any positive value of such 
that 

< ^ < ^, nl:= max{F(r, F)|r G (0, 1)}, 
m m 

there are two values of r, < vei < < 1 with Q"^ = mF{r, F). If F > 16/(9v^), 
there are no tetrahedral relative equilibria. This remark completes the proof. □ 



3. The characteristic polynomial 

The goal of this section is to obtain the characteristic polynomial, which will allow 
us to compute the spectrum of the Jacobian matrix corresponding to a tetrahedral 
relative equilibrium. The computations we perform and the conclusions we draw 
here will prepare the terrain for understanding the stability of the orbit relative to 
F and z. 

In equations ([2]), which describe the motion in rotating coordinates, the tetrahedral 
relative equilibrium becomes the fixed point 

(t) = (o)' (^2)^(0)' (|) = (v^/2)' (S)^(-vi/2)' 

e: = r/: = 0, ^ = 1,2,3,4. 

It is now convenient to introduce the linear operators Sd and So acting on 2 x 2 
matrices. 5*^ changes the signs of the elements on the diagonal, whereas So changes 
the signs of the other remaining elements. 



Stability in the positively curved 4-body problem 



Long but straightforward computations show that the Jacobian matrix correspond- 
ing to the vector field / corresponding to system (2) at the fixed point is given by 
the matrix 

■ / 
A B 



Df 



where 



5 = fidiag(5i,fi2,53,54), A = F + C, F = fi^ diag(Fi, F2, F3, F4), 



2 



Bo = 2 





-1 



B; = -r 



20;' 

/ mCii mCi2 mCi3 rnCii \ 



2 V -1-3^2 VS{l-z' 



-l+2z^ 




^3 4 V ^(1-4^2) _3 + i2^2 



, B4 = Sd{B3), 

, F4 = So{F3), 



C 



nil C21 XC22 XC23 XC24 
f^i C31 XC32 XC33 XC34 

y ?Tii C41 ^C42 ^C43 XC/n J 



12 ■ 



-2 
1 



7 C13- 



3^3 \ 



G 



21 ■ 



a 



4 V 3^3 -5 ) 

J( 

4Vy3(l + 2;22) i_g^2 
3 /-I + 9^2 _ ^g^4 ^3(1 _ ^2 



-2^2 
1 



_1/ 3-2^2 73(1 + 222- 

^3 + 

8 V 3y3(-l + 3^2) 5 _ 3;j2 



Ci4 — S'o(Ci3), 

, (^41 = 'S'o(C3i^ 



C24 — 5'o(C; 



23 J; 



3/ -924_g^2_^5 73(3;j4 + 4^2_i) 

''~8 V3V3(3^^-2z2 + i) -^z^' + Uz^-l 
/~i / 3 + 9^ 3^/3(32; 52; -|- 2)\ ^ , ^ \ ^ , ^ \ 

'-^34—^1 g 3 — 272"^ /' 42 — '-'o(,032j, <-/43 — Do(^034j 



_ 32 r< _ 3 



_3_ / -5 + 21^2 3^3(1-92 



16 V 3^3(1 - 922) 1 - 33^2 
So, we can write 

f]2^!!H ^3^^_3/2_ 
2 



1 - 1522 

-2 + 1222 J , 

2^ . 

C44 = 5'o(C33). 
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The eigenvalues of Df are the zeroes of the polynomial 

det(-C^/ + CB + A) = 0. 
Let us define fi such that ( = flfi. Then the characteristic equation becomes 

= det(-fiV^^ + f^^B + A) = 0. 

Let us introduce 

S := -n'^fi^I + + F + C = P(/i) + C, P(/i) := -fiV^^ + I^^B + F. 
Then 

P(/i) = fi2diag(Pi(/i),P2(/i),P3(/i),P4(/i)), 

where 



^3(/i) 



_ ^^(1 _ ^2) _ 1 + ^2 + ^2) + ^(1 _ ^2) 

-f (1 + 3Z^) + V^(^ - z') V ^ - - I + 3z^ 



and Pa{ij) follows from Pz{lj) by changing the sign in -\/3- Then p{jj) = det(S'). 

Before computing p{fi), it is convenient to perform some reduction and introduce 
additional notations. The first integrals associated to the energy and the SO (2) 
invariance give rise in to the factors /i^ and /i^ + 1, which we can ignore. To get 
further, recall first that, if the upper index denotes the transposed of a matrix, a 
2n X 2n matrix A is called infinitesimal symplectic if it satisfies the equation 

J A + A^J = 0, where ~ (^^j ^^'^ is the unit matrix. 

As the matrix P(yu) is infinitesimal symplectic, p{fi) contains only even powers of 
/i and, hence, we obtain with the notation M =: a simpler expression. We can 
further reduce the problem by considering a unique mass parameter. In general, 
we can discuss the stability in terms of the mass ratio F = mi/m, thus skipping 
the dependence on m. However, to study some limit cases, it will be also useful to 
consider e = m/mi instead of F. From now on we will simply denote the previous 
p{fi) by p{M), after changing the variable and skipping the factors M and M + 1. 

The characteristic polynomial p{M) has degree 6 in M, and its coefficients are 
polynomials of degree 8 in F that depend on z. The dependence on z is not of 
polynomial type due to the factors G~^^'^ and fi^. An important difference relative 
to the curved 3-body problem is that these factors cannot be easily "canceled" when 
multiplying by a power of Q"^, unless we take mi = 0. Introducing 

D = D(z) = G^'/^ = a(l + 3zY''^^ with a = 

9V3 
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the expression of p{M) becomes a huge polynomial, which can be fortunately sim- 
plified in part. 

Indeed, the factor F = AzVt^ appears in p{M) with multiplicity 3. Skipping it and 
further renaming the quotient as p{M), we obtain a polynomial of degree 6 in M 
whose coefficients have degrees 21 in z and 5 in D{z) and F. It is clear that the 
dependence on D can be decreased to degree 1, but then the degree in z increases. 
No other obvious factors appear. Whenever necessary, we will make the dependence 
on the other variables explicit by writing p{M, F, z, D{z)). 

As it is usually done in the 3-body problem, we can look for values of z and 
F related to bifurcations of the zeroes of p{M) that lead to changes in the spec- 
trum: either M = is a root or p{M) has a negative root with multiplicity at least 
equal to two. In the former case, after dividing by the factor z"^, the polynomial 
p{0,T, z, D{z)) has degrees 19, 5, and 3 relative to z,D{z), and F, respectively. In 
the latter case, after dividing by the factor F'^'^z'^^D{z)'^, the resultant of p{M) and 
-^p{M) produces a polynomial, denoted by Res(F, z, _D(z)), that has degrees 104, 
25, and 13 in z, D{z), and F, respectively. (Recall that if two polynomials P and Q 
have the roots ai, 02, . . . , and 61, 62, • • • , ^r;, respectively, then they have a common 
root if and only if Res(P, Q) = 0, where Res(P, Q) := nr=inj=i(^i ~ ^j) their 
resultant. In the present case M has to be seen as the variable of the polynomials 
and F and z as parameters.) Certainly, it can happen that p{M),p[0, F, z, D{z)), or 
Iies{r, z, D{z)) have some other non-trivial factor. But the dependence in D makes 
hard to recognize it. 

Hence, to study the stability problem, we will combine a numerical scan of the 
changes in the solutions Mi,i = 1, ... ,6, for some grids in F and z, with the the- 
oretical analysis done in the vicinity of some limit problems, which we will next 
introduce. 



Before proceeding with our numerical computations it is worth studying the be- 
haviour of the system in some simple limit cases, which we will later use to achieve 
our main goal of understanding the spectral stability of tetrahedral orbits. 

4.1. The restricted problem. If we take mi = 0, which is equivalent with F = 0, 
the matrix S has the block structure 



4. Three limit problems 




n^Piilj) +mC: 




11 




where 



11 



( 



n\i-n^) + ^z 



2/ifi2 



) 
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and Z(fi) is a 6 X 6 matrix such that all the terms have either a factor = 3mG 
or a factor X = mG~^^'^. Then 

det{S) = det(f}^Pi(/i) + mCu) det(Z(/i)). 

Note that from the matrix Z{fi) we recover the eigenvalues, and so the spectral 
stability of the Lagrangian orbits of the curved 3-body problem studied in [29] . 
These results, to be used in the next section, can be summarized as follows. The 
determinant of Z{fi) is a polynomial in M. After eliminating the factors M, + 1, 
and the exact solution given by Mq = —2^^(5—32;^) / (1+32;^), we obtain a polynomial 
of degree 3 in M, Q{M) (see also [22]), with polynomial coefficients in = 1 — 2^. 
In [2H] it was proved that there exist three values of r, < ri < r2 < < 1, where 
Hamiltonian-Hopf bifurcations occur, such that, for r G (ri,r2) U (r3, 1), the zeroes 
of Q are negative, and consequently those Lagrangian orbits for the curved 3-body 
problem are linearly (and not only spectrally) stable. For r e (0, ri) U (r2, r^), Q has 
a pair of complex zeroes, so the corresponding Lagrangian orbits are unstable. (For 
more details about Hamiltonian-Hopf bifurcations see [38].) 

In the restricted case, the stability of the zero-mass body located at (0, 0, 1) can 
be obtained by studying the matrix Q'^Pi{fi) + mCu. A simple computation shows 
that 

^ detin^Piifi) + mCu) = /i^ - c/i^ + (^^ + 2)', c = zG^^^ - 2. 

Then 

(5) M = /i^ = l(c ± v/-8(c + 2)) = i(c ± V-8^G3/2). 

If z > 0, /i^ becomes a complex number with real part different from zero. But if 
^ < 0, we obtain a couple of negative values for M, with an only exception that 
appears for c = -4, i.e. ^^(1 + 3^2^3 ^ 256/27 {z ^ -0.73176195875). For this z, one 
of the values of M is zero and the other value is negative. Of course, when z moves 
away from this exceptional value, the zero value of M becomes negative again. So, 
we can draw the following conclusion. 

Proposition 2. Considering the dynamics of the infinitesimal mass in the above 
restricted problem, the tetrahedral relative equilibrium is spectrally stable for negative 
values of z, but unstable for positive z. 

4.2. The 1 + 3 limit case: m = 0. As opposed to the previous restricted case, we 
now study the problem in which the body lying at (0,0, 1) is massive, whereas the 
other three bodies have zero mass. It is easy to see that 

det{S) = (fi2)8det(Pi(/i))det(P2(/i))det(P3(/i))det(P4(/i)) = (1^^) V^(l + /^^)^- 
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Skipping the trivial factors, the characteristic equation p{M) = reduces in the hmit 
to 

r(M) = M^(M + 1)^ = 0, 

which yields the roots /i = 0,/i = +i, and fi = — i, all of them of multiplicity 4. In 
other words, all the characteristic multipliers are equal to 1. 

This outcome is not unexpected. Indeed, if we use e = m/rrii as mass parameter, 
all the three equal masses are zero in the limit e = and their mutual influences 
vanish. Hence, the problem reduces to three copies of the 2-body problem, formed 
by the mass at the north pole and a body of zero mass. The changes in this highly 
degenerate situation for small e > will be studied in Section [61 

4.3. The solutions with z = 0. As we are also interested in the behaviour of orbits 
for small 2; > 0, it is also necessary to consider the solutions with 2 = 0. Skipping the 
trivial factors, we obtain again the limit equation T(M) = for all F. Again, this 
fact is not surprising because, for any positive F, we have that Q ^ 00 when 2 — 
and, therefore, the relative equilibrium requires larger and larger angular velocity. 
This means that the centrifugal force and the reaction of the constrains that keep 
the bodies on §^ are so large that the attraction of the mass lying at the north pole 
can be neglected. 

Regarding the cases in Subsections 14.21 and 14. 3^ we will further consider the be- 
haviour of the branches emerging from the solutions of T(M) = when the parame- 
ters e and z move away from zero. This analysis is cumbersome due to the presence 
of two parameters and of some long expressions. Furthermore, when e tends to zero, 
we want to study arbitrary values of z E (0, 1) and, when 2; — )■ 0, to consider arbitrary 
values of F in (0, 00). The bifurcations that occur in these cases will be studied in 
Section [61 

5. Numerical experiments 

The results of this section have been obtained using the polynomial p{fi) computed 
symbolically with PARI. According to the notation and reductions introduced above, 
we will also refer to this polynomial as p{M, F, z, D{z)). 

For given values of F and z, we first computed the zeroes Mi,...,M6 of the 
polynomial ]5(M, F, 2, D{z)). We used for the results plotted here a variable number 
of decimal digits, going up to 100 or more, and performed many additional checks. 

Recall that the complex zeroes, M, correspond to values of fi of the form ±a ± i/3, 
called complex saddles (CS); the real positive zeroes, giving values ±a for fi, are 
called real hyperbolic (H); and the negative zeroes, yielding ±i/3 for fi, are called 
elliptic (E). Changes in the stability properties occur when the zeroes pass from 
one type to another. The exceptional cases in which some zeroes of p{M, F, z, D{z)) 



14 



Florin Diacu, Regina Martinez, Ernesto Perez-Chavela, and Carles Simo 



are equal to zero or negative and coincident deserve attention to decide about the 
spectral stability of the solution, but they generically occur only in a zero-measure 
subset of (r, z). 

We will further use the coding E*H-'CS'^, where the exponents show the number of 
zeroes, M, of each type. Of course, the exponents satisfy the identity i + j + 2k = 6. 
In Figure [2] we display some numerical results. In the electronic version of this paper, 
the colour coding is 

E^^red, E^CS^^green, E^CS^^blue, E^H^ ^ magenta, E^H^CS^^pale blue. 

Hence, the observed transitions correspond to two types of bifurcations: 

- Hamiltonian-Hopf (for red green, green blue, and magenta — )■ pale blue) 
and 

- elliptic-hyperbolic (for red — magenta and green — >■ pale blue). 

The white zones are related to forbidden (F, z) domains, which correspond to < 0. 
In the printed version of this paper, the colours translate into grey shades as follows: 
red = black; blue = dark grey; pale blue = grey; magenta = light grey; green = very 
light grey. 

We first describe the case z > 0. The plot shows that for small F the orbit is 
unstable, except that near z = there is a line, emerging from F = 0, 2 = 0, where a 
super-critical Hamiltonian-Hopf bifurcation occurs and the system becomes totally 
elliptic. Again, for small F we find different zones for which one or two CS show up. 
As F — 0, the values of z at which the transitions occur tend to 

zi ^ 0.8299852976470169, Z2 ^ 0.7318602978602651, z^ ^ 0.3702483631504248, 

which correspond to the values ri,r2, and (see Section I^?T1) found in [29] for the 
Lagrangian orbits of the curved 3-body problem in S^. 

When F increases, as seen in the top left plot, a narrow red stable domain seems 
to persist near z = 0. The top right plot suggests that this is true up to F = 10. For 
some larger values of F, up to 10'^, this estimate seems to be still true. We can further 
ask about the limit behaviour when F — )■ cxd. The numerical evidence suggests, on 
one hand, that the limit value of z up to which the solution is totally elliptic is close 
to 0.03642; on the other hand, the boundary of one of the blue domains goes to z = 1 
and the domain disappears. The intermediate blue domain seems to shrink. Figure 
[3] provides more information: the blue domain shrinks to a point and increases again 
to merge with another blue domain born near F = 2.91, z = 0.822. It is remarkable 
that to the left of that point a tiny totally elliptic zone appears (one has to magnify 
the plot to see it). The blue domain for large F seems to tend to a limit width 
confined by values approaching 0.5 and ~ 0.94215. 

These numerical experiments raise the following theoretical questions for z > 0: 
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Figure 2. Regions of stability for the relative equilibrium orbits given 
as a function of the mass ratio F for the horizontal variable and z for the 
vertical variable (see the color code in the text). Top left: a general plot 
for F G (0,2]. Top right: magnification for F G (0, 10],^ G (0,0.12]. 
At the bottom we show a magnification of a narrow range F G (0, 0.1] 
for z < (left part) and an additional magnification around the little 
red triangle-shaped domain which can be seen near F = 0,2; = —0.8. 

(a) What happens when F — !■ oo (i.e. when e — )■ 0)? Do the red, green, and blue 
zones in Figure [3] on the right tend to a limit? 

(b) What happens for z very close to zero? In that case the value of f2 tends to 
oo when 2; — )■ and the limit is singular. A priori, some changes cannot be 
excluded in a tiny strip. 

(c) Which is the local behavior for F, z when both are positive and close to 0? 

We will return to these questions in Section O 
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Figure 3. Part of the domain z > for large values of T. The colour 
codes are described in the text. 

We further consider the case z < 0. As we already explained, the value of F is 
bounded by the condition = 0. In other terms, the boundary is parametrized by 
z G (—1, 0) as 

T{z) = T*{z) = -8^/[y3(l + 3z2)3/2]^ 

which we can also write as F*(z) = —9z{l + 3z^)-D/4. The most interesting domains 
appear for small values of F. There are two ranges of z, namely (— zi,— 22) and 
(— Z3, 0), in which the necessary conditions for linear stability are satisfied, in agree- 
ment with Proposition[2]and the corresponding results obtained in [29] for Lagrangian 
solutions. These ranges extend to small values of F; but there are two exceptions 
(both shown at the bottom of Figure [2]), namely when the axis F = is tangent to the 
red domains. These tangent points are located near z = —0.73892, z = —0.28396, 
and correspond to Hamiltonian-Hopf bifurcations. A red — )■ magenta transition is 
seen ending on a tangency to the vertical axis at (0, 0). The transition from the ma- 
genta to the pale blue domain seems also to be very close to the boundary F = T*{z) 
of existence of admissible values of z. 

These numerical results lead to the following problems in the case z < 0: 

(d) Prove that the transitions from stability to instability that occur in the re- 
stricted problem persist for F > 0. 

(e) Prove that there are exactly two additional values of z for which a curve of 
Hamiltonian-Hopf bifurcations is tangent to F = 0. 

(f) Analyze the vicinity of (F, z) = (0, 0) for z < 0. 

Like the questions (a), (b), and (c), we will address these problems in the next 
section. 
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6. The perturbation of the limit cases and the main result 

In this section we prove several results concerning perturbations of limit cases. The 
conclusions are summarized in Subsection 6.5. All proofs are analytical. The only 
use of some numerical information appears in the computation of the zeroes of a few 
polynomials of the form H{z, D{z)), a procedure that can be reduced to computing 
the zeroes of irreducible polynomials in z or by checking that some polynomials have 
a given sign at a given value of the variable. When we check that some polynomial 
is zero at a zero of some function, we either use the resultant or compute the zero 
with increasing number of digits. If d decimal digits are used and the zero is simple 
(respectively double), we check that the obtained value is zero up to approximately 
d (respectively d/2) digits. We increase up to a value that exceeds 1000. 

We begin with a lemma about the double zeroes of a function f{x,a,b), which 
depends nontrivially on two parameters a and b, i.e. neither nor fa nor are 
identically zero. In the applications to the present problem, x corresponds to the 
variable M, whereas a and b to z and F, respectively. We would like to see, for 
instance, if, for fixed F, two real negative zeroes of p that collide at a given value 
of z move away from the real axis, as well as what happens when F changes. The 
information we obtain is only based on the properties of /. We could exploit the fact 
that we are dealing with eigenvalues of an infinitesimal symplectic matrix (or a matrix 
conjugated to it), but some singular limit behaviour, such as when (T,z) — )■ (0,0), 
makes difficult to analyze perturbations of the limit case. Since we are interested in 
the vicinity of a point (x*, a*,b*), we shift the origin of the coordinate system to that 
point. We can now prove the following result. 

Lemma 1. Let f{x, a, b) be a real analytic function depending on the parameters a, b. 
Assume that for a = b = the function has a zero of exact multiplicity 2, located at 
X = 0, i.e. /(O, 0, 0) = /x(0, 0, 0) = and, for concreteness, /xx(0, 0, 0) > 0. We want 
to study the behaviour of f in a neighbourhood of (0, 0, 0). For fixed b = 0, we have: 

(i) ///a(0,0,0) > when a increases, crossing the value a = 0, the roots move 
away from the real axis. The case /a(0, 0, 0) < is similar when a decreases. 

(ii) ///a(0, 0,0) = 0, consider faa{0, 0,0) and fxa{0, 0,0). If the discriminant 
fxa ~ fxxfaa O't (0, 0, 0) is positivc, thc roots remain real. 

Let now b vary. Then: 

(a) Under the assumptions of (i), there exists a line a = h{b) along which f 
has double zeroes in the x variable, and when a increases, crossing the value 
a = h{b), the roots move outside the real axis. 

(b) Under the assumptions of (ii), if fh{0, 0,0) > 0, there exists a curve b = k{a), 
with positive quadratic tangency to b = at a = 0, such that the zeroes of f 
pass from real to complex when crossing the line b = k{a). 
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(c) Under the assumptions o/ (ii), and z//f,(0,0,0) = 0, there are two curves, 
say hi{b),h2{b) (eventually complex or coincident), tending to {a,b) = (0,0) 
when b ^ 0. If they are real and distinct, say hiib) < h2{b), then the zeroes 
of f are real if a < hi{b) or a > h2{b) and complex if a & {hi{b), h2{b)). 

Proof. The cases (i) and (ii) are elementary, since the Newton polygon in x, a involves 
the vertices (2, 0) - (0, 1) and (2, 0) - (0, 2), respectively. 

To prove (a) we can assume that fxx{0, 0, 0) = 1, /a(0, 0, 0) = 1, scale the variables, 
and write f{x, a,b) = + a + 0{x^, ax, a^) + bg{x, a, b). To find a double zero, we 
can use the Implicit Function Theorem to express x as a function x = x{a, b) from 
the equation fx{x, a, b) = 0. Inserting this x in the equation f{x, a, b) = 0, we obtain 
a relation between a and b. The Implicit Function Theorem allows us then to express 
a as a function of b. 

To prove (b), we scale the variables and apply a linear change in the (x, a) variables, 
after which we can write that f{x, a, b) = x"^ — + ct)!^) + ^(1 + ^(^) o, b)). 

From the equation /j. = we obtain x = x{a,b) as in (a), and if we insert x in /, 
we can write 6 as a function of a that starts with a positive quadratic term in a. 
Undoing the linear change and scalings simply deforms the picture linearly. 

To prove (c), we proceed as before, obtain x = x{a, b) and insert it in /. But now 
the linear term in b is absent, while there is a nonzero quadratic term in a. The 
existence of the two branches follows by using a Newton polygon in a, b. □ 

Remark 1. In the exceptional case of item (c) in which hi{b) = h2{b), the function 
/ can be written, after an eventual shift of x, as f{x, a, b) = (x^ — h{a, by)f{x, a, b) 
with /(O, 0, 0) 7^ 0, and h{a, 6) = if the parameters a and b satisfy a = hi{b). Then, 
for values a,b with h{a,b) = 0, x has a double zero and the relations f{0,a,b) = 
0, fx{0, a, b) = 0, fa{0, a, b) = 0, and fb{0, a,b) = hold. 

In particular, the resultant of f = 0, fx = with respect to x, gives h = with 
multiplicity 2. However, the resultant of fx = 0, fb = gives h = with multiplicity 1. 
Note that in the case (a) the resultant of /x = 0, = is far from zero along h = 0. 
In the case (c) with hi{b), h2{b) real and distinct, the resultant of fx = 0, fb = gives 
a single line with 6 as a function of a with multiplicity 1. 

6.1. Analysis of the case F — t- oo. In this case only z > has sense. Let us 
consider p{M) in terms of e = 1/F. The term in £° is M^{M + 1)^, whereas the 
terms in ,j > 0, are polynomials of degree 6 in M that, in turn, have polynomials 
in z,D{z) as coefficients. 

To discuss how the double zero M = bifurcates as a function of e, we compute 
the Newton polygon in the e, M variables. After simplifying by a numerical factor, 
we obtain 

16M2+eM((1728^^-720^^ + 722)D-192^^+4822)+e2[(466562^°-38880^^ + 11988^^ 
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-1620z^+81z'^)D^+{-10368z^+6912z'^-1512z^+108z^)D+576z^-288z^+36z^]=0, 

a quadratic equation in M with discriminant zero. Hence the dominant term of M 
is of the form 

Mo{e) = ^z{4z^ - l)h{z,D{z)), h{z, D) = {3 - 18z^)D + 2z. 

It is easy to check that the factor h in Mo{e) is positive in the interval z G (0, 1). 
Indeed, h is positive at 2 = and at 2 = 1 and dh/dz has only two zeroes in the 
interval (0, 1) at which h is positive. Therefore the only value of z at which Mo(£) 
changes sign in the interval (0, 1) is at 2 = 2:5^0 = 1/2- Hence, Mo(e) becomes positive 
for 2 > 1/2, giving rise to instability, in agreement with the lower bound of the blue 
domain for large F obtained in Section [51 However, the analysis up to now shows that 
the two branches emerging from M = are real and coincide. It could happen that 
higher order terms take them away from the real axis even for values of z G (0, 1/2). 

As usual, we introduce another variable by the transformation M = MQ{e)+eN. 
Substitution into the characteristic equation and division by e"^ gives the dominant 
terms in the new Newton polygon. They turn out to be, up to a numerical factor, 
of the form 

A^^ + e3h%z'^ - l/4:)g{z, D), g{z, D) = h{z, D)k{z, Df, k{z, D) = D{9z^ - Qz) - 1. 

We note that g{z, D) is positive. Therefore 

M{e) = Mo{e) ± e"^/^ ^/3H^{1/ A - z^)g{z,D). 

This shows that the roots M evolving from zero are real, negative, and distinct for 
z < 1/2, and are complex with positive real part for z > 1/2. When the variable z 
crosses the value 2; = 1/2, a Hamiltonian-Hopf bifurcation occurs. 

Let us analyze the solutions evolving from the quadruple solution M = —1. We 
introduce a new variable, which we denote again by A^, such that now M = —l+eN. 
The lower order term in e is a term in whose coefficient Q{N, z) is a polynomial 
of degree 4 in A^ with polynomials in z, D{z) as coefficients. 

First we would like to determine the behaviour of the function for z small. The 
dominant terms are 

8A^^ - IS^aA^^ + lOS^^aA^^ + 1701z^a^N - 78732z^a^ 

where, we recall, a = D{0). A Newton polygon method tells us that the dominant 
terms in the solutions are 

A^i = 9za/4, A'2 = 62^ N3 = -36z^a, A'4 = 81z^a/4. 

In particular all the solutions are simple and negative, ensuring local spectral stability 
near z = 0. 
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To study the behaviour of the function for larger values of 2, we compute the 
resultant of Q{N,z) and -^Q{N, z). After skipping some powers of z and D, we 
have a polynomial of degree 40 in z and degree 10 in D. It is easy to check that this 
polynomial has only two simple zeroes for z G (0, 1), located at 

^4,0 ~ 0.036420258329089021, ze,o ^ 0.942152758989663983. 

At 

2:4^0 a couple of roots meet and become complex, whereas at zq q these roots return 
to the real domain. 

Hence, we can now summarize the bifurcations for e small as follows. 

Proposition 3. Fore = l/F tending to zero, there exist three functions, z^le), z^{e), 
and ZqIe), tending to 2:4,0,^5,0? '^i^'d zq^, respectively, at which Hamiltonian-Hopf bi- 
furcations occur. They are sub-critical at z^le) and super- critical at z^le) and ZQ{e). 
Therefore the character of the fixed points is of type for z G (0, 24(e)), of type 
E^CS^ for z G {z,{e), z,{e)) U {ze{e), 1), and of type E^CS^ for z G {z,{e), ze{e)). 

Proof. The above analysis and items (i) and (a) of Lemma [1] complete the proof. □ 

6.2. The case of small positive z. To study this case it is convenient to replace 
D{z) by a(l + 32^)^^/^, as done before, and to expand the binomial up to the required 
order. First we study the solutions emerging from M = 0. 

We proceed as in the previous subsection, by regarding M as a function of z for z 
around 0. First we find two branches that coincide at order 1 in 2;: M = —daez/A, 
where we use again e to denote F"^. Then we seek the terms in z"^ that are also 
coincident. At the third step there appear two branches in 2;^/^ with opposite signs. 
Summarizing, the solutions evolving from M = are 

M{z) = -9aez/4: + {Sla'^e^ - 24e)zyi6 ± 9{ae^y^^ z^/y2. 

That is, the two values Mi,M2 emerging from M = are real negative and they 
only differ in the 0{z^^'^) terms. For further reference we denote them by Mi (with 
+) and M2 (with -). 

For the solutions that evolve from M = — 1, we write M{z) = — 1 + A^ and compute 
the Newton polygon in the z, N variables. After simplifying constants, the dominant 
terms in the polygon are 

4A^^ - 9azN^e + bAaz^N^e^ + 35842^iVe^ - 165888a2^e^ 

From this expression we obtain the dominant terms of the four branches, already 
separated at this first step, 

Ni = 9aez/4:, N2 = Gez^, N3 = 256e^V(3a), N4 = -mGez^ {27 a). 

We can now summarize the above results as follows. 
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Proposition 4. For any fixed value of the mass ratio T in (0, oo), there is a range 
of values of z close to zero, the upper limit of the range depending on T, such that 
the spectral stability of the orbit is preserved. 

This result explains the red domain displayed in the previous figures for z > 
small. Note, however, that, as soon as F — )■ 0, we have e — > oo, and then the range 
of validity of Proposition S] is not uniform, since it can go to zero when F — )■ 0. 
This fact has been already put into the evidence in the Figures at the corner near 
{T,z) = (0, 0) in the first quadrant, where a bifurcation line is seen to emerge from 
(0,0). The required analysis follows in next subsection. 

6.3. Study of the vicinity of {T,z) = (0,0) for both z > and z < 0. When 
approaching (0, 0) in the (F, 2)-plane, we have a singular problem. Depending on the 
direction, the value of can tend to any real non-negative value. Therefore, before 
proceeding with the analysis, we must add a short description of the difficulties we 
face, based on the following numerical experiment. 

For this purpose, we wanted F to be neither too small (to exclude a configuration 
with too close lines), nor too large (to be inside the domain of interest), and thus 
chose F = 0.03. Then the values of the solutions M were computed as a function of 
z. The value of z corresponding to fi^ = is z* ^ -0.00649642410717306. Figure 
m on the left plots three of the M values, after multiplying them by which is 
a measure of the distance from z to z*. For the values which appear to be almost 
constant and tend to coincide when z — )■ {z*)~ , we also changed the sign. This means 
that when approaching = 0, one value of M seems to tend to +00, whereas two 
values seem to tend to — cx). 




-0.015 -0.012 -0.009 -0.006 -0.015 -0.012 -0.009 -0.006 0.005 0.01 0.015 0.02 0.025 



Figure 4. Evolution of the solutions Mj as functions of z for F = 0.03 
as an illustration of the behaviour of the zeroes of p{M, F, z, D{z)) for 
small z and F. See the text for a detailed explanation. 

The other three values are plotted in real scale on the middle plot. We can clearly 
see the Hamiltonian-Hopf bifurcation near z = —0.00793. After that value of z, we 
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only plot the real part of these solutions. For a previous value near z = —0.01072 
there is also a double root, which seems to avoid a bifurcation. Furthermore, a 
multi-precision study in narrow intervals provides evidence that, when z comes very 
close to z*, the solutions that became complex have a real part that turns to positive 
(without giving rise to a bifurcation), whereas the negative one remains negative, 
tending to a value, when z ^ z*, that tends to zero when F — 0. In conclusion, we 
can expect two double zeroes in that domain, only one giving rise to a bifurcation. 

When we pass to values oi z > 0, the behaviour of the zeroes is shown in Figure 
m on the right. The solution Mi evolving from remains real while M2 provides 
a Hopf bifurcation when it collides with the solution that starts at — 1 + A^2- The 
branches starting as — 1+A^3 and — 1+A^4 remain real and stay close to —1. Similarly, 
the solution starting as — l + A^^i remains real, although it coincides with the solutions 
Ml and M2 at two values of z, which are very close. Hence, we should expect three 
double roots for that value of F at some small z values, only one giving rise to a 
bifurcation. 

The above numerical evidence will help us obtain some theoretical results. But 
before proceeding with the analysis of the bifurcations, it is useful to study the 
behaviour of the function in the vicinity of the curve = near {T,z) = (0,0). 
This approach presents interest in itself because it gives us the eigenvalues and allows 
us to interpret the results we will later obtain. 

For this purpose, we follow an approach different from the one that led us to Figure 
m Instead of fixing F and allowing z to vary, we fix 2; < near and vary F. Recall 
that T*{z) has been defined as the limit value corresponding to Q"^ = 0. We further 
define 7 = F*(2;) — F and want to study what happens when 7—7-0'^. Our findings 
are summarized below. 

Proposition 5. For z < close to zero, let T*{z) be the value ofT corresponding to 
= and 7 = F* — F. Then, when 7 — O"*", the roots of the characteristic 
polynomial, for a fixed value of z, and then when 2 — 0~, behave as follows: 

(i) Two roots are real and negative, and when multiplied by ^/\z\ they tend to a 
common value Xi{^)> behaving close to the limit, when 7 — t- 0''", like Xii^) =t 
0(^/72). The value of Xiiz) tends to — 8/-\/3 when 2 — 0^ with a dominant 
term that is linear in z. 

(ii) A third root is real and positive, and when multiplied by ^/\z\ it tends to a 
value X3(^); from below, behaving close to the limit like X3(^) — ^{.iz)- The 
value of Xsi^) tends to 8/\/3 when 2 — 0^ with a dominant term that is 
quadratic in z. 

(iii) A fourth root is real negative and tends to a value X4(^) linearly in 7. The 
limit value X'^i.^^ tends to —1 when z— )-0~ with a dominant term linear in z. 
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(iv) The last two roots are complex, and when multiplied by ip = 'y / z"^ they tend 
to a common nonzero real value. The real part multiplied by ip tends to 
Xbi^); which tends linearly in z to 27/2 when 2 -> 0~. The imaginary parts 
multiplied by y/^p tend to values ixei^); which tend linearly in z to ±-\/54 
when z — )■ 0~ . 

Proof. We only sketch the main steps, the full result following then easily. The 
characteristic polynomial p is written as a function of z (a polynomial in z and 
D{z)), 7, and M, and we look at the Newton polygon in the variables 2;, 7, both of 
them to be seen as small, involving the exponents (0, 5), (1, 4), (2, 3), (3, 2), (7, 1) and 
(11,0), with coefficients that are polynomials in M, given by 

--f^2^°M^{M + 1)^ + -f^za2^3'^M{M + lf{M + 2) + -f^{zaf2'^3^{M - 1)(M + 1)^- 

-f\zaf2^3%M + 1)3 + 7zV2^3^(M - 1)(M + 1) - z^^a^2''3^\M + 1). 

The last equation has the obvious solution M = — 1, which is simple. By adding the 
remaining terms in p, the solution mentioned in (iii), denoted as M4, is obtained. 

Using the side between (0,5) and (3,2), the variables z and 7 should be of the 
same order. This suggests to change the variable M to by M = Nza/'y. Setting 
7 = 0, and simplifying by some powers of z and A^, we obtain the polynomial 
64Ar3 _ ;l44jy2 _ 324^ + 729^ ^hich has the roots N = 9/4 (double) and N = -9/4. 
Hence, we obtain solutions Mi^2 whose main terms are 9az/{4'y), as stated in (i), 
and M3 whose main term is —9az/{A'y), as stated in (ii). 

Using the side between (3,2) and (11,0), we obtain that 7 = 0{z^), which suggests 
a change of variable from M to defined as M = Nz'^/'-f. As before, setting 7 = 0, 
and simplifying the powers of z, we obtain the polynomial 4A^^ — 108A^ + 729, which 
has the roots 27/2 (double). Hence, we obtain solutions Ms g whose main terms are 
27^^/(27), as stated in (iv). 

Let us denote by Mf^\i = 1,...,6, the first approximations of the 6 solutions 
found up to now. As usual, we write Mi = Mf^'' + AM^*-"'' and substitute them in the 
initial equation. We then obtain the new Newton polygons and find the corrections, 
as described in the statement. □ 

Remark 2. The properties described in Proposition [5] agree with the observed fact 
that the points close, but below, the line fi^ = belong to the pale blue domain. 

We return now to study a vicinity of (F, z) = (0, 0). We know that the bifurcations 
we are looking for are associated to M = or to double roots. We begin with the 
case M = 0. 

By skipping a suitable factor, the polynomial p(0, F, z, D{z)) has 
32T^ + lUDzT^ + 162D^z'^T + 3645^^2^ 



24 



Florin Diacu, Regina Martinez, Ernesto Perez-Chavela, and Carles Simo 



as Newton polygon in the z,r variables, where its coefficients still depend on D{z). 

The last two terms give the branch whose dominant term is F = which 
can be written around ^ = as r = — 45a2;'^/2 = — 802;^/a/3, in perfect agreement with 
the numerical results. This phenomenon is easily identified with the red to magenta 
transition in Figure |2] on the left, both top and bottom. The other vertices give rise 
to a factor T = -9zD{l + 3z^)/A + 27z^/2, double up to order 4 in z, but which 
is located between z = and the curve = and, hence, outside the admissible 
domain. 

We will further study the double roots. As mentioned at the end of Section [31 the 
resultant polynomial Res(r, z, D{z)) is huge, but for (2;, V) near (0, 0) it is still feasible 
to compute the Newton polygon. The relevant vertices bounding the three sides of 
the polygon have exponents (16,0), (12,2), (9,4), and (0, 13). After simplifications, 
the first side from (16, 0) to (12, 2) leads to branches with dominant terms given 
by the solutions of the equation 3242;^a^ + = 0. They are complex and can be 
discarded. 

The second side, with endpoints (12, 2), (9, 4), gives the condition for the dominant 
terms bAz^a — V"^ = 0, with the real solution z = (4x3^/^)~^F^/'^. It is easily identified 
as the curve that separates the green and red domains from each other in Figure [2] 
top, both left and right, near (0,0). 

The third side gives branches with dominant terms of the form F = 9za/4,r = 
— 9za/4, and F = — 36za/25 of multiplicities 4, 3, and 2, respectively. We begin 
with the case of multiplicity 3. Setting F = z{—9D{z)/4 + 7) in the resultant and 
simplifying by constants and powers of z and D, we obtain the polynomial 

1287^ + 3888D^z-f'^ + 52488D^z^-f + 177U7D^z^. 

This gives raise to one branch which, to order 2 in z, and using the full D{z) function, 
not onlyD(O) = a, is of the form F = —z9-D(z)/4—243z^-D^/8, easily identified as the 
transition from magenta to pale blue. The other root is double, 7 = — 27 z'^D/ 4. At 
the next step, F = z{— 9 D{z)/ 4)+ z^{— 27 D/ 4+ ip), we obtain again a double solution, 
•0 = 27z/8. But the first part of F given by —9zD{z){l + 3z^)/4 is the boundary 
= 0. Hence, the obtained double branch is already outside the admissible domain 
at order 4 in z. 

We now consider the branch beginning with F = — 36za/25 of multiplicity 2. In 
fact, the successive Newton polygons that we computed in the expression of F as a 
power series in z always give multiplicity 2. Hopefully this branch of double zeroes 
corresponds to the double zeroes that appear in Figure HJ in the middle, and do not 
give rise to a bifurcation. We further computed two additional resultants. Up to 
now we are using Iies{r, z, D{z)), obtained from the elimination of M between p 
and dp/dM. After simplification, the degrees in z,D{z), and F are 104, 25, and 13, 
respectively, as mentioned before, and the polynomial contains 6779 terms. 
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Let Res2(T, z, D{z)) be the resultant from p and dp/dT and Res^iT, z, D{z)) the 
resultant from dp/dM and dp/dT. The corresponding degrees and numbers of terms 
are similar (111, 27, 11, and 7453 for Res2 and 96, 23, 11, and 5474 for Ress). But 
the important thing is that the Newton polygons of Res and Res2 give, up to the 
computed order, a branch of multiplicity 2, while the one of Resa is simple. More 
precisely, the Newton polygon of Ress has degree 11 and factorizes as 

{Ay - 9f{Ay + 9)^(8?/ - 9){25y - 36)(25y + 36)(80y2 + 720y + 1377), 

where y denotes the ratio r/{zD). We are interested in the ratio y = —35/25, simple 
as claimed. A few terms of the expansion of F as a function of z, D[z) are obtained, 
in a recurrent way, as 

^ , 36L>, ^3-(1944L)2 + 54) ^4 + (186624L)^+31968L'2_^144) 

^ = "^^^-^)' = ' 625 ' = ' 25D ' 

c/5-(35831808D6 + 59222259L'^-12906L)2^768) 



50D 

^6 + (3439853568L'S+10318220802L)6^42691698D^+226512D2^2048) 



9a 

25D 

Hence, the branch is double and, according to Lemma [Hand Remark [T], no bifur- 
cation occurs along that line. As a side information we note that along that double 
branch, for z < small, the value of M is close to —1. 

Finally we consider the branch starting with F = 92;a/4 of multiplicity 4. Writing 
F = 2;(9-D/4 + and substituting in Res, we obtain the polynomial 

-16^^ + 288^^^ - 13682^^^ + 6482^^2 - 81^^ 

which factorizes as (4(7| — 362;5'2 + 92^)^. Hence, the terms of order 2 in z give rise to 
two double solutions, with coefficients 9/2 ± 3\/2. As in the previous case, keeping 
with Res, we obtain double solutions for the computed next terms. 

We will further use Ress. As mentioned before, one of its factors is 4?/ — 9 with 
multiplicity 3. Setting F = z{9D/4+g2), we obtain the Newton polygon in the {g2, z) 
variables as 

-512^^ + (8748L>^ + 6656)z^^ + (-78732^^ - 19584)^^^2 + (19683L>^ + 4608)z^ 

which factorizes as {Agl - 3Qzg2 + 9z'^){-128g2 + (2187^2 + 512)^). The last factor 
is irrelevant for our purposes and the quadratic factor gives the two branches with 
dominant terms g2 = {9/2±3\/2)z, which are simple. Hence, as in the previous case, 
the two branches of Res are double and they give rise to no bifurcation. Additionally, 
we can mention that these double solutions occur for M close to —1/2 and that from 
the plot in Figure H] on the right we expect them to be close. 

We can now summarize the results obtained in this follows. 
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Proposition 6. In a vicinity o/F = 0, 2 = 0, for F > there are three lines giving 
rise to bifurcations, all emerging from (0,0).' 

(i) A line of E ^ H transition, for z <{], having a cubic tangency with the axis 
F = 0. 

(ii) A line of E"^ — j- CS transition, also for z < 0, which has a quadratic tangency 
with the line corresponding to = 0. 

(ill) A line of E"^ — t- CS transition, for z > 0, for which z is of order F^^^. 

6.4. Analysis of F near zero. In this subsection, we need only to consider bifur- 
cations that occur away from a vicinity of (F, z) = (0, 0), since the behaviour in the 
neighbourhood of this point has been already studied in the previous subsection. 

For positive z, we must only show that the changes of stability that occur for 
the 3-body problem persist when we add the small mass mi. As already shown in 
Proposition [21 the behaviour of the small mass gives instability for the full 4-body 
problem. Using the results in [29] and of Lemma [U^a), it follows that the changes of 
stability of the 3-body problem persist in the case F > small. 

Next we pass to the more interesting case z < 0. The changes of stability found in 
the curved 3-body problem also persist when we pass to F > 0, according to Lemma 
[Uj^a), and the stability of the body of mass zero does not change the stable domains 
for F = 0. However, new changes can occur when passing from F = to F > if 
some of the additional zeroes of the form ([5]) for the restricted problem coincide with 
one of the curved 3-body problem. 

Figure [5] shows all the relevant zeroes, when real, simultaneously as function of 
z. The exact solution Mq (see Section Hj) can be identified as the curve starting at 
(—1, 1) and ending at (0,0) (see also formula (7) in [22]). Double zeroes involving 
Mq should not be taken into consideration in the light of the explanations given in 
that paper. In Figure we identified two double zeroes that are the responsible for 
the Hamiltonian-Hopf bifurcations observed in Figure [2l 

To locate them, we consider the polynomial Q{M) introduced in Section 14. 1[ It 
is convenient to express the coefficients of Q{M) in terms oi Z = z"^ = 1 — r"^. We 
will further denote this polynomial by Q{M, Z). Also we write the equation for the 
additional zeroes ([5]) as S{M) = 0, where 

S{M)= iM^ + 2M+—^ + lj -z^G%M-lf. 

Using G = 3(1 -|- 3z'^)/A, we can write S{M) as a polynomial in M with polynomial 
coefficients in Z to be denoted as S{M, Z). Then we compute the resultant TZ{Z) of 
Q{M, Z) and S{M, Z) to eliminate M. The polynomial TZ{Z), of degree 29, factorizes 
as TZ{Z) = 7li{Z)7l2{Z) , with factors of degrees 14 and 15, respectively. A part of 
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Figure 5. A plot of the negative values of /i^ as a function of z < 
for the limit problem F = 0. The large dots show the location of the 
double zeroes associated with the Hamiltonian-Hopf bifurcations that 
can be seen as tangencies of the red domains with F = shown in 
Figure |2] bottom, left and right. The curves of zeroes associated with 
the restricted problem are easily identified as having a tangency with 
z = at n^ = —l. The branch reaching the lower boundary of the plot 
continues down up to ^ = — 1. 

the expressions is 

7^l(Z) = 31381059609Z^^ + 135984591639Z^3 + . . . + 843283683Z - 2985984, 

n2{Z) = 1162261467Z^^ + 5036466357Z^^ + . . . + 133996544Z + 16777216. 

The polynomial TZi{Z) has exactly 4 real zeroes with Z G (0,1), approximately 
located at the following values of z: 

-0.071519103755, -0.114735617843, -0.330240264422, -0.736842605000, 

whereas 7^2 (^) has exactly two zeroes in the same interval corresponding to the 
values of z 

Zi ^ -0.7389177458229170, Z5 ^ -0.2839588732787964. 

This fact is in agreement with the plots in Figure To show that the two zeroes 24, z^ 
for F = give rise to a bifurcation, it is enough to check that c}Res(F, z, D{z))/dT ^ 
at the points (0, Zj, D{zj)),j = 4, 5. The computed values are ~ —0.626660386126 at 
Z4 and ~ 27.6667376231 at z^, far away from zero. 

In the case of the other four double zeroes, we obtain values of (9Res(F, z, D{z))/ dT 
equal to zero (with the expected accuracy, see the beginning of the present section). 
Imposing the condition of double zero for M, dp{M,r, z, D{z))/dM = 0, and sub- 
stituting it in p{M,r, z, D{z)), we obtain, to low order, a double branch of double 
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zeroes in the {T,z) variables. Using now Resa, as we did in Subsection I6.3[ we ob- 
tain that the branch is simple. There is no need to employ the Newton polygon; 
the Implicit Function Theorem is enough because the linear coefficients are nonzero. 
Hence, no bifurcation related to the zeroes of TZi occurs. 
We have thus proved the following result. 

Proposition 7. In the passage from the restricted to the general problem for z < 
and a small mass ratio T, excluding a neighbourhood o/(r, z) = (0, 0), already studied 
in Proposition changes in the stability properties occur along lines of the (F, z) 
plane. Three of these lines tend to the values —zi, ~Z2, —z^ when F — j- 0. Additional 
changes occur along two curves, with quadratic tangencies to the line T = at the two 
values Z4 and z^ where the characteristic polynomials of the curved 3-body problem 
and the restricted problem have zeroes in common. In all these cases a Hamiltonian- 
Hopf bifurcation occurs. 

6.5. The main result. We can summarize the conclusions obtained in the above 
propositions by saying that they validate the numerical results near the relevant 
boundaries of the domain (F,^;), i.e. near F = 0, F = oo,z = 0. We have found all 
bifurcations produced by perturbation of the limit problems. These properties also 
describe a general view on the problem of stability of tetrahedral orbits in the curved 
4-body problem in S^. The properties rigorously proved above are now summarized 
by the following result. 

Theorem 1. We consider the tetrahedral solutions of the positively curved A-body 
problem on with a fixed body of mass mi located at the north pole and the other 
three bodies of equal mass m located at the vertices of an equilateral triangle orthog- 
onal to the z-axis. Let be T = — > 0. Then: 

m 

(1) For r — )■ oo there are three functions, 2;4(1/F), 2;5(1/F), and zq^I/T), tend- 
ing to z^^Q, z^fi, and zq^, respectively, at which Hamiltonian-Hopf bifurcations 
occur. The respective tetrahedral relative equilibria are 

- spectrally stable of type for z G (0, Zi{l/V)). 

- unstable of type E^CS^ for z e (z4(l/F), 25(1/F)) U {zq{1/T), 1), and of 
type E'CS' for ze{z,{l/r),z,{l/r)). 

(2) For any fixed value of the mass ratio F in (0, oo), there is a range of values of 
z of the form (0, 5(F)), i(F) > 0, such that in that range the orbit is spectrally 
stable. The value of z{r) tends to 24,0 for F — 00 and 2(F) — )■ for F ^ 0. 

(3) For any z > fixed when F approaches the orbit is unstable. 

(4) For any z < fixed and close to zero, let T*{z) be the value corresponding to 
Q"^ = and 7 = F* — F. Then, when 7 — )■ 0^, the corresponding tetrahedral 
relative equilibria are unstable. 
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(5) Around the point T = 0, z = there are six sectors, ai,a2,cr-s,(T4,a5,aQ, or- 
dered counterclockwise, in which the type of the orbits are: E^, E^H-^, E^H^CS^, 
no solutions, E^, and E^CS^ , respectively. The dominant terms in the bound- 
aries of the sectors are of the form F = 0, F = —ciz^, z = — a/3F/8 — C2F^, 
z = — a/3F/8 — csF"^, z = 0, z = C4F^/^, and F = 0, the first four with z < 
and last two with z > 0. All the coefficients Ci are positive. 

(6) For small F and z < there exist five curves, ipiy4'2,4'3,4''i, and ip^, at 
which stability changes of Hamiltonian-Hopf type occur. The first three are 
transversal to the line F = 0, whereas the other two are tangent. These curves 
reach F = for the following values of z: 

zi ^ -0.829985, Z2 ^ -0.731860, Z3 ^ -0.370248, 

Zi ^ -0.738918, Z5 ^ -0.283959. 

For small T, in particular, the tetrahedral relative equilibria are spectrally 
stable for z between: ipi and the lower branch of ip^; the upper branch of ip^ 
and ip2; "^3 and the lower branch of ipb,' and between the upper branch of 
and the curve bounding a\ in item (5), above. 

7. Conclusions and outlook 

In this last section we will draw some final conclusions about the stability of 
tetrahedral relative equilibria and propose three problems that, in order to be solved, 
would require certain refinements of the methods we applied here. 

We can now summarize the stability results we obtained in this paper about the 
relative equilibria of the tetrahedral 4-body problem in by displaying the full 
bifurcation diagram. To complete the above analysis of the limit cases, we present 
the diagram computed from the resolvent Res(F, 2, and from the conditions 

p(0, F, 2;, = 0. In both cases, given a value of z^ we obtain a polynomial 

equation for F. We computed the zeroes numerically and discarded the ones that do 
not give rise to any bifurcation. We checked the facts that occur here by looking at 
the derivatives with respect to M and z at the solutions found. Figure E] depicts the 
results. As horizontal variable we used F/(l + F) in order to display the full range 
of F G [0,cx)]. 

A possible continuation of the present work is the study of the linear stability for 
pyramidal solutions, i.e. orbits of the positively curved n-body problem, for n > 4, 
with a fixed body of mass m\ located at the north pole and the other n — 1 bodies 
of equal mass m lying at the vertices of a rotating regular polygon, orthogonal to 
the 2;-axis. But the analytic methods pursued here have limits. It seems that the 
symbolic computations and the related analysis would become insurmountable for 
n larger than 7 or 8. Even a purely numerical study must be done very carefully. 
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Figure 6. Bifurcation diagram for the relative equilibrium solutions. 
Variables displayed: (r/(l + T),z). With the only exception of the 
boundary of the admissible domain = for z < and the line 
going from (0, 0) to the tip of the = line, which corresponds to 
E-H bifurcation, all the other lines correspond to Hamiltonian-Hopf 
bifurcations. 

Another interesting problem is to analyze the linear stability of tetrahedral orbits 
in Finally, the stability of tetrahedral orbits in §^ when the ^-coordinate of the 
three equal masses is not constant, but varies periodically in time, would also be a 
problem worth approaching. 
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